'''
Created on Oct 3, 2010

@author: oabalbin
'''
import numpy as np
import scipy as sp
from scipy import stats
#from collections import deque, defaultdict
import RNAseq.io.translate_gene_names as tn
from RNAseq.common.classes import myArray

def get_gene_anotations():
    '''
    returns ucsc annotation based on knowGene, kgXref, and knownIsofoms table.
    change so you do not harcode the paths
    '''
    inputfile_genes = '/data/ucsc_tables/knownGene.txt'
    inputfile_Hugo = '/data/ucsc_tables/kgXref.txt'
    inputfile_isoforms = '/data/ucsc_tables/knownIsoforms.txt'
    knownGenes = tn.read_known_genes(inputfile_genes)
    knownHugoNames = tn.read_hugo2ucsc_names(inputfile_Hugo)
    knownIsoforms = tn.read_known_isoforms(inputfile_isoforms)
    
    return knownGenes, knownHugoNames, knownIsoforms


def print_expression_matrix_stdfiltered(thisGeneArray, outfile_matrix, lstd, gct_format=False, gene_dictionary=[]):
    """
    """
    
    ucsc_kg, hugo_kg, ucsc_known_isoforms = get_gene_anotations()
    if not gene_dictionary:
        print "Translating ucsc to hugo names and ref seqnames"
        #gene_dictionary = tn.get_isoform_list(ucsc_known_isoforms,thisGeneArray.gendict, hugo_kg)
        gene_dictionary = tn.get_hugo_names(thisGeneArray.gendict, hugo_kg)
 
    
    list_of_samples = thisGeneArray.get_list_of_samples()
    #print list_of_samples
    list_of_genes = thisGeneArray.get_list_of_genes()
    list_of_intervals = thisGeneArray.genInterval
    indicators = np.array(range(len(list_of_genes)))
    
    gene_std = thisGeneArray.sdv_array(thisGeneArray.expVal)
    gene_med = thisGeneArray.median_array(thisGeneArray.expVal)
    std_perc = sp.stats.scoreatpercentile(gene_std,lstd)
    med_perc = sp.stats.scoreatpercentile(gene_med,lstd)
    
    ngenes_excluded = len( set(indicators[gene_std < std_perc]).union(set(indicators[gene_med < med_perc])) ) 
    ngenes_included = len(list_of_genes) - ngenes_excluded
    print ngenes_excluded, len(gene_std[gene_std < std_perc]), len(gene_med[gene_med < med_perc])
    print ngenes_included
    if gct_format:
        outfile_matrix.write("#1.2\n")
        outfile_matrix.write(str(ngenes_included)+'\t'+str(len(list_of_samples))+'\n')
        outfile_matrix.write("Name"+'\t'+"Description"+"\t"+",".join(list_of_samples).replace(',','\t')+'\n')
    
    else:
        outfile_matrix.write("Name"+'\t'+"Description"+"\t"+",".join(list_of_samples).replace(',','\t')+'\n')
        
    #for i in range(thisGeneArray.expVal.shape[0]):
    for g, i in thisGeneArray.gendict.iteritems():
    
        if (gene_std[i] < std_perc) or (gene_med[i] < med_perc):
            continue
        eline = list(thisGeneArray.expVal[i,:])
        # especially for proteomics experimens        
        nzeros = eline.count(0.0)
        #if nzeros >= np.ceil(0.70*len(eline)):
        #    continue
        #eline = list(thisGeneArray.expVal_foldchange_mean[i,:])
        chr, loc = list_of_intervals[g].split(':')
        glength = int(loc.split('-')[1]) - int(loc.split('-')[0])
        outfile_matrix.write(gene_dictionary[i]+'&'+str(i)+'&'+str(list_of_intervals[g])+'\t'+str(glength)+'\t'+",".join(map(str,eline)).replace(',','\t')+'\n')
        #outfile_matrix.write(list_of_genes[i]+'_'+str(i)+'\t'+str(glength)+'\t'+",".join(map(str,eline)).replace(',','\t')+'\n')


def print_expression_matrix_isoforms_stdfiltered(thisGeneArray, outfile_matrix, lstd, gct_format=False, gene_dictionary=[]):
    """
    write the expression matrix summarizing the isoforms
    """
    
    ucsc_kg, hugo_kg, ucsc_known_isoforms = get_gene_anotations()
    if not gene_dictionary:
        print "making the isoform dict"
        gene_dictionary = tn.get_isoform_rows(ucsc_known_isoforms,thisGeneArray.gendict, hugo_kg)
    
    list_of_samples = thisGeneArray.get_list_of_samples()
    list_of_genes = thisGeneArray.average_isoforms(gene_dictionary)
    indicators = np.array(range(len(list_of_genes)))
    list_of_intervals = thisGeneArray.genInterval
    
    gene_std = thisGeneArray.sdv_array(thisGeneArray.expVal_isoaverage)
    gene_med = thisGeneArray.median_array(thisGeneArray.expVal_isoaverage)
    std_perc = sp.stats.scoreatpercentile(gene_std,lstd)
    med_perc = sp.stats.scoreatpercentile(gene_med,lstd)
    
    print len(list_of_genes), len(gene_std)
    
    ngenes_excluded = len( set(indicators[gene_std < std_perc]).union(set(indicators[gene_med < med_perc])) ) 
    ngenes_included = len(list_of_genes) - ngenes_excluded

    if gct_format:
        outfile_matrix.write("#1.2\n")
        outfile_matrix.write(str(ngenes_included)+'\t'+str(len(list_of_samples))+'\n')
        outfile_matrix.write("Name"+'\t'+"Description"+"\t"+",".join(list_of_samples).replace(',','\t')+'\n')
    
    else:
        outfile_matrix.write("Name"+'\t'+"Description"+"\t"+",".join(list_of_samples).replace(',','\t')+'\n')
    
    for i in range(thisGeneArray.expVal_isoaverage.shape[0]):
        
        if (gene_std[i] < std_perc) or (gene_med[i] < med_perc):
            continue
         
        eline = list(thisGeneArray.expVal_isoaverage[i,:])
        #nzeros = eline.count(0.0)
        #if nzeros > 0.4*len(eline):
        #    continue
        
        outfile_matrix.write(list_of_genes[i]+'\t'+str(i)+'\t'+",".join(map(str,eline)).replace(',','\t')+'\n')

def print_expression_matrix(thisGeneArray, outfile_matrix):
    """
    """
    list_of_samples = thisGeneArray.get_list_of_samples()
    #print list_of_samples
    list_of_genes = thisGeneArray.get_list_of_genes()
    outfile_matrix.write("geneName"+'\t'+",".join(list_of_samples).replace(',','\t')+'\n')
    
    for i in range(thisGeneArray.expVal.shape[0]):
        eline = list(thisGeneArray.expVal[i,:])
        outfile_matrix.write(list_of_genes[i]+'&'+str(i)+'\t'+",".join(map(str,eline)).replace(',','\t')+'\n')

def print_foldchange_matrix(thisGeneArray, outfile_outliers_matrix):
    """
    """
    list_of_samples = thisGeneArray.get_list_of_samples()
    #print list_of_samples
    list_of_genes = thisGeneArray.get_list_of_genes()
    outfile_outliers_matrix.write("geneName"+'\t'+",".join(list_of_samples).replace(',','\t')+'\n')
    
    for i in range(thisGeneArray.expVal_foldchange.shape[0]):
        eline = list(thisGeneArray.expVal_foldchange[i,:])
        outfile_outliers_matrix.write(list_of_genes[i]+'\t'+",".join(map(str,eline)).replace(',','\t')+'\n') 
